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- - - We study optimization of superconducting properties of inhomogeneous nanowires. The main goal 

' of this research is to find an optimized geometry that allows one to maximize the desired property 

I of superconductors, such as the maximum value of local superconducting gap or total condensation 

energy. We consider axially symmetric design of multi-layered nanowires with possibility to adjust 
and change the layers thickness. We use numerical solution of the Bogoliubov-de Gennes equations 
to obtain the local superconducting gap for different arrangements of the inhomogeneous structures. 
, The value of the optimized properties can be up to 300% greater compared to a non-optimized 

geometry. The optimized configuration of multilayers strongly depends on the desired property one 
wants to optimize and on the number of layers in the nanowire. 
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I. INTRODUCTION 



o, 

O , Recent advances in fabrication and experimental measurement of spatially inhomogeneous superconducting ma- 
' terials open up new opportunities in optimal design of targeted properties of superconducting materials. Optimal 
Oh' design is an approach where one calculates the properties of a given structure or device and then optimizes the shape, 
^ I composition or mutual orientation of the parts of the structure with the purpose of achieving some function, some 
target property of the structure. With proper chosen geometry (that determines boundary conditions for quantum 
"j^ i systems) one can expect enhancement of the target property at certain location, or even observe a "quantum mirage" 
in confined corral-like nanoscale systems^. 
I , The ideas of optimal design are following the same line of thinking as the band structure engineering^!^, or prediction 
• of the complex materials structure using information-theoretical methods-. Notion of optimal design also can be 
Ch viewed as an extension of similar ideas to the design of decoherence-protected quantum state engineering in quantum 
^ computing^ 

I I [ In this work we propose to extend the notion of design of quantum states to the case of superconducting materials. 
We will focus on the superconducting nanostructures. Examples of the specific parameters that one might want to 
optimize in superconducting materials are local superconducting gap A(r) and maximum of the condensate energy 
^ Econd C!C / c?r|A(r)p. These are the main target parameters we will focus on below. Specifically, we are going to 
investigate the interplay between strong spatial inhomogeneities, confinement and geometry effects in superconducting 
nanoscaled systems. To elucidate this interplay we will focus on the influence of the spatial arrangement of layers of 
different superconducting materials on the target parameters we mentioned above. 

Obviously, the search for optimal design can only be performed using physically correct and computationally acces- 
sible model. In our approach we use the Bogoliubov-de Gennes equations which is well-suited to describe inhomoge- 
neous superconductors microscopically. Since three-dimensional spatial inhomogeneity is not tractable without using 
, a supercomputer, we resort on the case of cylindrical quasi-one-dimensional systems, considering axially symmetric 
L| ■ nanowires. The axially symmetric superconducting wires were recently fabricated and have shown their usefulness for 
. ^ ' prospective measuring devices^i^iiSiiiS. 

Before we proceed with the results we will make a few technical comments. Spatial inhomogeneity destroys trans- 
Jh lational symmetry and can create localized states on a certain length scale. If the inhomogeneities are created by 
. . 1 disorder, it usually may be characterized by simple, low order, and spatially localized correlation functions. However, 
this approach cannot be applied in the case of engineered, ordered inhomogeneities. In this case inhomogeneities are 
characterized by spatially delocalized correlation functions of high order. As a result, this problem is not tractable 
analytically. One needs to use numerical solutions with a fine mesh and at the same time allow for variations in the 
solution as we search for optimized configuration. This approach was not possible computationally ten years ago and 
only now becomes feasible. 

In our approach we abandon the usual assumption that the designed inhomogeneity can be treated as a small 
perturbation. Detailed analysis can demonstrate that the size of the space of accessible solutions grows with the 
maximum allowed amplitude of the inhomogeneities. Thus, a relatively big inhomogeneity may bring new effects 
and functionalities, which are not accessible by small perturbations. The assumption of relatively large spatial 
inhomogeneities place our search for optimal configurations beyond the realm of applicability of the Anderson's 
theorem. The Anderson's theorem for non-magnetic alloys and zero external field states that to a first approximation 
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the excitation spectra are the same for alloy as for the pure metalii. This theorem resorts on the assumptions of 
small spatial inhomogeneity of the pair potential A(r) and chemical similarity of the introduced impurities. Both of 
the assumptions do not hold in configurations that turned out to be optimal. 

If designed inhomogeneities are large, a naive strategy to maximize the local superconducting gap A(r) may be 
achieved by controlling quasiparticle local density of states at and near the Fermi level. On a deeper level one may 
expect the interplay between correlations of quasiparticles on the scale of the coherence length ( and the scale of 
spatial inhomogeneity. Note that the coherence length may be also modified by the spatial inhomogeneity, thus the 
problem of the finding of the proper length scales should be solved self-consistently. 

Important effect that comes into play is the confinement effeclji^, when a quasiparticle subband appearing due to 
the finite size quantization happens to be close to the Fermi surface. In this case the density of states is increased, 
that leads to enhancement of the pairing potential. 

The pairing potential A(r) may also depend of the finite size of the system via the phonon quantization. However, as 
it was shown for thin film&12,, this effect is a minor correction to the quasiparticle quantization mechanism mentioned 
above. Here we leave consideration of the phonon quantization effects for future publications and consider phonons 
to be the same as in the bulk material. 

This paper is organized as follows. In Section II we present an approach, based on numerical solutions of the 
Bogoliubov-de Gennes equations. In Section HI we discuss the results of numerical solutions and present main 
results. In Section IV we discuss main qualitative finding of the optimal design approach and its applicability to other 
problems. We conclude in Section V with the discussion of main results and the outlook of future approaches. 

II. APPROACH AND BASIC FACTS 

In this work we use Bogoliubov-de Genes formalism^**. It is a mean- field formulation of microscopic theory of 
weak coupling superconductors, and it is particularly well suited for spatially inhomogeneous problems. It is also 
appropriate for the relatively small hcterostructures in the clean limit that we are interested in. 

Using the microscopic theory for inhomogeneous superconductor, one can explore the rich structure in the local 
density of states (LDOS), which now can be measured directly by using STMi^. Such STM studies, unlike the previous 
experimental studies^^ which explored only spatially averaged properties, may provide valuable information about 
local variations of the pair potential A(r). 

We are going to use real space discretization techniqueii, that unlike using basis of Bessel functions for the cylindri- 
cally symmetric problem a^^'^^ , leads to the eigenproblem with sparse matrices. In this case the eigenproblem solution 
can be found much faster using sparse iterative solvers. In this work we choose the wire to be hollow to simplify 
numerical calculations and to avoid uneven real space descretization for the numerical treatment of the singularity 
near the origin. However, our results will not be changed qualitatively in the case of solid wires of the same size. 

The Bogoliubov - de Gennes (BdG) equations for the quasiparticle amplitudes m„ and w„ with excitation energy 
En > can be written in a compact form 

- . . 

(-TT^V^ - ^i)a,^Pn{r) + A(r)a^Vn(i-) = -B„Vn(r), (1) 

where we denote V'n = {un,Vn)'^, and ax,d'z are the Pauli matrices. Here /x is the chemical potential. Any effects of 
the band structure of the material are included in the quasiparticle effective mass m*. In the presence of magnetic 
field one should replace the Laplace operator by (p — f A)^. 

The superconducting pair potential A(r) for s-type superconductor must be determined self-consistently from the 
solutions of the BdG equations, as: 

A(r) =5(r)^u„(r)<(r)[l-2/(£;„)]e(fit^z,-i?„), (2) 

n 

where f{En) is the Fermi distribution function, lojj is the Debye cut-off frequency, Q{) is the step function, and g(r) 
is the inhomogeneous coupling constant. The summation is performed over the indices n corresponding to £"„ > 0. 
Throughout this work, we use chemical potential /i = and temperature T = OK. The coupling constant is chosen in 
order to obtain the bulk value A^uik ~ 10 K. Note that as in other studiesi^, similar model parameters were chosen 
for illustrative purposes only, and were not intended to reproduce realistic materials. 

Assuming z axis to be an axis of the axial symmetry in the wire of length L and imposing periodic boundary 
conditions along this axis, the solution of Eq. ([T|) has the form: 

7/'„(p, 0, z) = d„(p) exp(?;0)/\/27rexp(ifczz)/VL. (3) 
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The BdG equations for the radial part can be written as: 

(-^— - ^i)&zan{p) + A{p)a^an{p) = EnUnip), (4) 

where we use the notation 

vl = dydp^ + p-^d/dp-r'p-''-k^, (5) 

for the radial part of the Laplace operator. Indices I and — 27rn^/L are the azimuthal quantum number and the 
wave vector in the z direction correspondingly. 

To compute the self-consistent solution to both Eqs. (|l|2p . we assume an initial guess Ao{p) — A^uik for the order 
parameter. Then by means of numerical diagonalization, we compute the solutions Ui{p), Vi{p), for each combination 
{I, kz}, corresponding to Aq{p). We use these solutions to generate new guess for the order parameter. 

In order to avoid instabilities during the self-consistent calculations, that we found typical for strongly inhomoge- 
neous geometries, we use a simple mixing procedure: 

A™„(p) = aA,(p) + (l-a)A,_i(p), (6) 

where a is an adjustable parameter, < a < 1. Initially we set a = 0.5. To insure convergence, we increase the current 
value of a by 5% if the relative deviation between two consequent steps S"* = max^ |Ai(p) — A,i_i(p)|/ maxp |Ai(p)| 
has decreased, S'- < And we decrease a by 20% in the opposite case, > S'-~^. 

The computed A„ij^(p) is used for the next iteration step. 

We repeat iterations until we achieve the acceptable level of accuracy (S*' < 10^'^). After the end of the procedure, 
we perform an additional step with a = 1 to ensure convergence of the obtained solution. It usually takes 50 — 100 
iterations to converge. 



III. RESULTS 



A. Homogeneous wires 

In order to check our simulations we first reproduce the bulk limit by considering a relatively big system. We 
assume a hollow homogeneous cylindrical wire with zero value boundary conditions on the inner and outer surfaces. 
For such small systems the role of the boundary conditions, as we see later, is crucial. Zero boundary conditions we 
chose does not imply necessarily contact with a normal metal. It can be vacuum (or just air) outside the wire. 

We choose the outer radius of the wire Rmax — 50 nm, and the inner radius as Rmin = 0.5 nm (see Fig. [T]). The 
length of the wire is set to be L = 500 nm with periodic boundary conditions that is a sufficient approximation for 
an infinitely long wire. The temperature in all our simulations is set to zero T = K. 

We find that the radial part of the self-consisted pairing potential A(p) is close to the bulk value with some deviations 
near the boundaries. The slight rise of the pairing potential can be attributed to the local approximation for the 
pairing potential in the BdG equations, so it has a nature of the Gibb's phenomenon. Now let us compare a hollow 
homogeneous cylindrical wire with inner radius Rmin and outer radius Rmax, assuming zero boundary conditions, and 
a homogeneous film of thickness D = Rmax — Rmin- The results of the calculations are shown in Fig. [21 The spatially 
dependent superconducting order parameter is normalized to the maximum value Aq for the slab. We chose Rmax—^^ 
nm and i?mm=l-5 nm. For all our simulations we use the length of the wire L = 500 nm with periodic boundary 
conditions. The maximum value of the A(p) for the hollow wire is about 30% larger than for a slab. Clearly, for the 
wire the superconducting order parameter has asymmetric shape due to different local curvature at different values of 
p. Next we calculate A{p) for wires of different thickness. In Fig[3]we plot the superconducting order parameter A(p) 
for the fixed outer radius Rmax = ^5 nm and variable inner radius Rmin — ^-^, 3.75, 7.5 nm. We found that max{A(p)} 
monotonously decreases as a function of Rmim together with decreasing of the local curvature. 



B. Inhomogeneous wires 

In this section we focus on few geometries of inhomogeneous nanowires. In the first geometry we consider a hollow 
axially symmetric nanowire composed of two layers of two different superconducting materials. We model different 
materials by using different coupling constant g. We have chosen the thicknesses of the layers to be equal, i.e. the 
first arrangement is g{p) = gi, for Rmm < R< {Rmm + Rmax)l^, and g{p) = g2, for {Rmin + Rmax)/^ < R< Rmax- 
The second arrangement one can obtain by switching the constants gi and g2- 
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Now let us demonstrate an unexpected effect of the geometry on the superconducting order parameter of an 
inhomogeneous system. We consider a wire which has an inner coaxial layer of a better superconductor with the 
coupling constant gi — 2g and outer layer with the coupling constant 92 — 9- For the second arrangement we 
consider the inverse sequence of the layers, namely, the wire has inner coaxial layer of superconductor material with 
the coupling constant gi = g and outer layer with the coupling constant 52 — 25. The first choice is depicted in Fig[4] 
by solid line, and second choice - by dashed line. Note that the calculated A(p) is normalized to maximum value Aq 
for the slab geometry of the same thickness D = Rmax — Rmin- Note, that because of the spacial variation of the 
local curvature the exchange of the superconducting layers leads to nonequivalent results. Interestingly, although for 
the second arrangement the system has a bigger volume of the better superconductor, than in the first arrangement, 
the maximum value of the pair potential A{p) in this case is smaller. This effect is due to local enhancement of the 
quasiparticles coupling in the effective potential U^ff oc for p — > 0, originated from the angular momentum term 
in the radial part of the Laplacian Eq.®. 

Now let us consider an inhomogeneous wire composed of three axially symmetric layers. In Fig. [5] we compared 
three different arrangements for the three layer system. We have chosen the layers to be of the same thickness, i.e. 

dip) = ffl: for Rmin < R < (Rmax + 2i?™„)/3, g{p) = 52, for {Rmax + 2i?mi„)/3 < R < (2Rmax + Rmin)/S and 

g{p) = (73, for (2Rmax + Rmin)/-^ < R < Rmax- The coupliug coustaut for the one of the layers (either gi or 52 or 
53) is set to 2g, and the two others are set to g. For this type of inhomogeneity the best choice that maximizes the 
superconducting order parameter A(p) is the one, where the better superconductor is in the middle. Since the layers 
considered here are thinner compared to the layers shown in Fig. [H the superconducting order parameter of the layer 
closest to the inner surface is significantly reduced by the proximity effect of the nearest boundary with A(i?mm) = 0. 

Now let us consider an optimal design problem. It is known, that the condensation energy E^ond is an important 
quantity that characterizes superconducting systems. Obviously, for inhomogeneous system it may vary for different 
arrangement of the superconducting layers. It is interesting to find an optimal configuration that maximizes Econd- 

In order to isolate volume effects we assume a superconducting layer of a fixed volume, that is equal to one third of 
the total volume of the wire: V2g = Tr{Rmax ~ ^min)^/'^- This layer has the couphng constant 2g. The rest material 
in the wire has the coupling constant g. The problem is to find the best placement of the layer to maximize the 
condensation energy E^ond 

nRma. 

Econdoc / A^{p)pdp. (7) 

Let the position of the layer's boundary, that is closest to the axis of the symmetry labelled by i?, then the layer's 
boundary closest to the outer surface of the wire will be Rr = •\/V2g/ (tt-L) — R"^ where L is the length of the wire. 

In Fig. [6] we plot the condensation energy Econd normalized to its maximum value Emax , as a function of its 
position R. For our simulations we choose Rmin — 1-5 nm, Rmax = 15 nm, L = 500 nm. One can clearly observe that 
the condensation energy has a maximum near i? « 5 nm. For this placement the reducing effect from the boundary 
A(i?mm) = is minimized, while the effective enhancement of the superconducting order parameter due to the local 
curvature is still taking place. There are some other local maxima of the condensation energy, that we attribute to 
finite size resonance^^. 

As we mentioned above, it is now possible to measure the local density of states (LDOS) directly using STM 
technique. In the STM experiments^ LDOS is proportional to the local differential tunneling conductance. In Fig. 
[7] we plot the calculated local density of states po for a homogeneous superconducting nanowire that is determined 
through 

Po{p,E) ^Y.^\ur,{p)\H{E - En) + \vn{p)?5{E + E^)]. (8) 

n 

We again consider a hollow wire with Rmin = 1-5 nm and Rmax = 15 nm. The length of the wire is chosen to be 
500 nm. We use E = Abuik for these calculations. In Fig. [7] one can see relatively big spatial variations of the local 
density of states due to the confinement of quasiparticles in the cylindrically symmetric boundaries. One possibly 
may observe these variations by scanning the end face of the wire by a STM tip. 

IV. DISCUSSION 

The notion of design of quantum properties is well established. The examples of this quantum state optimization 
include band structure engineering^'-, engineering of complex materials structure's. Notion of optimal design also can 
be viewed as similar to the ideas of decoherence-protected quantum state engineering in quantum computing^i^. 
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We proposed to apply these ideas to the superconductors as one example of macroscopic quantum state. The 
optimal design approach allows us to investigate the optimal structures that maximize the functional property of 
nanowires, like maximal condensate energy and superconducting gap. 

We find that depending on the desired property, specific design is different. Specifically wc find that for two layered 
superconducting wires maximal gap is obtained when the superconductor with the stronger pairing potential is placed 
on the inside. The reason for this enhancement is a confinement effect in the cylindrical geometry. 

In the case one wants to optimize the maximum critical current one would need to place stronger SC on the outside 
because this will be the region where the screening current fiows the most. We explicitly demonstrated that depending 
on the function one wants to optimize the superconducting wire design with the same ingredients can be very different. 
On the other hand, in the case of three layer superconductors we see that the largest gap (with enhancement factor 
of about 300%) is obtained when the strongest superconductor is placed in the middle of a sandwich. 

Optimal design ideas presented here are applicable to other configurations. One might look at the multilayered 
structures obtained by the molecular beam epitaxy and other methods. Here the combination of superconducting and 
insulating or magnetic layers (ferromagnetic, multiferroic and other layers) would allow us to investigate the optimal 
design of multilayered structures with competing interactions. One particularly interesting case is the experiment by 
Bosovic et aP^. In their multilayered configuration pseudogaped normal regions are alternating with the supercon- 
ducting regions in multilayers. The unusual proximity effect at very large distances of up to 100 Ahas been observed. 
Wc can take these configurations as a guidance on how possible future multilayered structures might look. The ideas 
of optimizing the competing interactions in these configurations are becoming relevant. 

V. CONCLUSION 

We consider spatially inhomogeneous superconducting materials of the axial symmetry. We also studied the inter- 
play between local curvature and spatial inhomogeneity of the system. We find that due to the cylindrical geometry, 
different arrangement of superconducting layers are not equivalent because of the different local curvature, and the 
interchange of the sequence of the layers gives different results. We find that placing of the better superconducting 
material layer closer to the symmetry axis results in the higher values of the pairing potential. However, it is not 
true in the case when this layer is thin enough. In this case the boundary effects will reduce effect of the small local 
curvature. We have proved this for the case of three layers. 

We also studied the optimal design problem of placing of the better superconducting layer to achieve maximum 
pairing potential. The optimal placement is a result of compromise between the local curvature effects and effect of 
the close boundaries. For more complicated geometries and larger number of layers one should perform search with 
the help of global optimization techniques, like Genetic Algorithms. The obtained solution arc most likely will not be 
accessible with any known analytical methods. 

For ultrathin nanostructures we considered in this work, localized impurity effects and phase fluctuations may 
play a significant role. The impurities and defects can be modeled on the mean field level by introducing localized 
impurity potentials, similar how it was done, for example. The phase fluctuations effects may compete with the 
pairing potential enhancement due to quantum confinement of quasiparticles in the cylindrical geometry of a nanowire. 
However, such effect goes beyond mean fleld BdG theory used in this work. The effect of quantum fluctuations may 
be an interesting topic to study in the future. 

As a possible way to extend our research, we consider to study nonequilibrium transport properties of inhomogeneous 
nanowires (super-current distribution, etc.). Since the nanowires under our considerations are much smaller than the 
characteristic field screening length, the inhomogeneity of the pairing potential will not play a significant role for 
relatively small fields and currents. However, for larger currents, closer to the critical value, the geometry and size 
effects may significantly change nonequilibrium characteristics of inhomogeneous superconducting nanowires, resulting 
in higher values of the critical current. It may also be interesting to study how the critical current value depends on 
the length of nanowire, and how the pairing potential changes along the nanowire of a finite length. 

We believe, the optimal design has a bright future in nanoscience since it can potentially convert the available 
computational power into improvement of target properties of rationally structured materials and may substantially 
increase efficiency of engineered devices. 
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APPENDIX 



We already outlined in the introduction the main points of our approach. Here we will elaborate on a more technical 
details of our numerical approach. In our simulations the radial part of the Bogoliubov-de Gennes equations Eq. ^ has 
a singularity at the origin p = that needs uneven discreization mesh and mixed type of the boundary conditions'^. 
In order to avoid these complications, we restricted ourselves to consider hollow cylinders with inner radius Rmin and 
outer radius Rmax- We assume Dirichlet type of the boundary conditions, namely: 

Then we introduce a transformation 

a{p) = b{p)/^p (A-2) 
that removes in Eq.Q the term with the first derivative with respect to p: 

("T^^W - {l^ " l/4)p-' - kl - p)aXip) + Aip)&Xip) - eX{p). (A-3) 
2m* 

The transformation Eq. (jA-2[l guarantees that the descretized Hamiltonian in Eq. (jA-3[) is represented by a symmetric 
matrix, that ensures much faster numerical diagonalization. 

The radial part of the Bogoliubov-de Gennes equations Eq. (jA-3|) is discretized on a equidistant mesh in real 
space. For this purpose we use the 13th order 12-point discretization rule for the second derivative^. The exact 
diagonalization of the resulting sparse matrix is then performed numerically using the ARPACK eigensolver— available 
freely on the web^3,. We performed a comparison of the numerical results and the analytical solution available for 
A(p) = 0, and find a good agreement using the number of discretization points N = 200. 

For every iteration step of the self-consistent calculations, the numerical diagonalization procedure is performed for 
different values of the wave vector and the azimuthal quantum number I. The procedure stops when no eigenvalues 
lay in the Debye window determined by huo- 
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FIG. 1: The geometry of the system is depicted in Fig. la and is represented by an axially symmetric hollow wire. The calculated 
radial part of the superconducting order parameter A(p) (see Fig. lb), is normalized to the bulk value, for a uniform infinite hollow 
superconducting wire with inner radius Rmin = 0.5 nm and outer radius Rmax = 50 nm. The length of the wire is chosen to be 500 nm. 
For our calculations we assume temperature T = K. 
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FIG. 2: The geometries of the systems under consideration are depicted in Fig. 2a and are represented by a hollow wire and a thin 
film. The calculated radial part of the superconducting order parameter A{p) (see Fig. 2b), for a uniform hollow superconducting wire 
with inner radius Rmin = 1-5 nm and outer radius Rmax = 15 nm (solid line). The length of the wire is chosen to be 500 nm, and 
the superconducting order parameter for uniform superconducting film with the thickness D = 13.5 nm (dashed line). The two other 
dimensions of the film are set to 500 nm. For our calculations we assume temperature T = K. 
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FIG. 3: The calculated radial part of the superconducting order parameter A(p), for uniform infinite hollow superconducting wire with 
the inner radii _Rmi„=1.5, 3.75, 7.5 nm and the constant outer radius Rmax = 15 nm. The length of the wire is chosen to be 500 nm Note 
decrease of the order parameter with the decreasing of the thickness. For our calculations we assume temperature T = K. 




FIG. 4: The geometry of the system is depicted in Fig. 3a and consists of 2 axially symmetric layers. The order parameter A(p)(see 
Fig. 3b) for a nonuniform hollow wire with Rmin = 1.5nm and Rmax = 15 nm. The length of the wire is chosen to be 500 nm. For our 
calculations we assume temperature T = K. The thick solid line corresponds to the case when the layer with stronger coupling constant 
2g closer to the center, dashed line corresponds to the case when the layer with stronger coupling constant close to the outer surface. The 
weaker superconductor has the coupling constant g. Each layer thickness 6.75 nm. Note the non equivalence of the arrangements due to 
the local curvature. 
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FIG. 5: The geometry of the system is depicted in Fig. 3a and consists of 3 axially symmetric layers. The order parameter A(p) (see Fig. 
4b) for a nonuniform hollow wire with three layers and three different arrangements, the better superconductor has a doubled coupling 
constant. Rmin = 1.5nm and Rmax = 15 nm. The length of the wire is chosen to be 500 nm. For our calculations we assume temperature 
T = K. 




FIG. 6: The geometry of the system under consideration is similar to what is shown in Fig 5a. The better superconducting layer is in the 
middle. The condensation energy E oc f A?{p)pdp of a heterostructure as a function of the position of the better superconducting layer. 
The better superconductor has a doubled coupling constant 2g, the worse superconducting layers have g. Rmin = i--5nm and Rmax = 15 
nm. The length of the wire is chosen to be 500 nm. For our calculations we assume temperature T = K. 
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